What you did since the beginning of the year (with link functions or not) was something like this \[
y = X\beta + \epsilon
\] Where \(y\) is the response variable, \(X\) the set of predictors and \(\epsilon\) the error term.
These models, - assume that all variables are directly observed/manifest - allow measurement errors only in endogenous variables - are just particular cases of SEM
SEM formula
In fact, the structural model of a SEM (i.e., excluding latent variables) can be expressed with the following equation: \[
Y = X^\ast B' + \zeta
\] Where - \(Y\) is the (n x p) matrix of endogenous variables - \(X^\ast\) is the n x (p + q) matrix of endogenous and exogenous variables - \(B\) is the (p + q) x (p + q) coefficient matrices - \(\zeta\) is the (p x q) matrix of errors in the equations This looks pretty similar to the regression formula, but with some matrices! Univariate regression models are just a special case of this formula where the parameter matrix is full of 0!
Imagine you want to predict scores in the test we will do at the end of this corse (\(y\)), based on your prior statistical knowledge (\(x_1\)) and interest (\(x_2\)): - First define the model
Second simulate the data
A first example with simulated data
# Simulate knowledge and interest as predictors of Yset.seed(12)N =100x1 =rnorm(N)x2 =rnorm(N)y = .35*x1 + .20*x2 +rnorm(N)d <-data.frame(x1,x2,y)
cor(d)
x1 x2 y
x1 1.0000 0.0159 0.386
x2 0.0159 1.0000 0.118
y 0.3863 0.1185 1.000
#cov(d)
x1 x2 y
x1 0.7484 0.0138 0.364
x2 0.0138 0.9982 0.129
y 0.3644 0.1291 1.189
lm regression
# fit a regression modelm <-lm(y ~ x1 + x2, data = d)summary(m)
Call:
lm(formula = y ~ x1 + x2, data = d)
Residuals:
Min 1Q Median 3Q Max
-2.8022 -0.6244 -0.0259 0.7150 1.8090
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0251 0.1009 -0.25 0.80
x1 0.4846 0.1172 4.14 7.5e-05 ***
x2 0.1226 0.1015 1.21 0.23
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.01 on 97 degrees of freedom
Multiple R-squared: 0.162, Adjusted R-squared: 0.145
F-statistic: 9.36 on 2 and 97 DF, p-value: 0.000191
Introducing lavaan
lavaan (latent variable analysis) is actually THE package for SEM. You can use it to estimate a wide family of latent variable models, including: factor analysis, structural equation, longitudinal, multilevel, latent class, item respons, and missing data models…
..But also simple regressions
Model fit and info
library(lavaan)ml <-"y ~ 1 + x1 + x2"#1 + gives the interceptfit <-sem(ml, data = d)# summary(fit, rsquare=T)
lavaan 0.6-19 ended normally after 1 iteration
Estimator ML
Optimization method NLMINB
Number of model parameters 4
Number of observations 100
Model Test User Model:
Test statistic 0.000
Degrees of freedom 0
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
[...]
# And we can plot itlibrary(semPlot)semPaths(fit, whatLabels ="parameters",edge.label.cex =1.5, rotation =2,residuals = F,sizeMan =10,curve =1.9,edge.color="black", edge.label.color="black")
Model matrices
#We can also look at the matrices#The parameters matrixinspect(fit)$beta
y x1 x2
y 0 2 3
x1 0 0 0
x2 0 0 0
inspect(fit, "estimates")$beta
y x1 x2
y 0 0.485 0.123
x1 0 0.000 0.000
x2 0 0.000 0.000
#The residual var-covar matrixinspect(fit)$psi
y x1 x2
y 4
x1 0 0
x2 0 0 0
inspect(fit, "estimates")$psi
y x1 x2
y 0.987
x1 0.000 0.741
x2 0.000 0.014 0.988
Basic lavaan syntax
As you can see, the regression syntax of lavaan is actually the same as lm, but there is much more in lavaan.
Model specification sintax:
Syntax
Function
Example
~
Regress onto
Regress B onto A: B ~ A
~~
Residual (co)variance
Variance of A: A ~~ A Variance of A and B: A ~~ B
=~
Define a reflective LV
F1 is defined by items 1-4: F1 =~ i1 + i2 + i3 + i4
<~
Define a formative LV
F1 is defined by items 1-4: F1 <~ i1 + i2 + i3 + i4
:=
Define non-model parameters
u2 := x + y
*
Label or fix parameter
Z ~ b*X labels the regression as b
Basic lavaan functions
Function
Command
sem() / cfa()
Fit the SEM model (cfa is nested in sem…which is nested in lavaan)
fitMeasures()
Return fit indices of the SEM model
inspect()
Inspect/extract information that is stored in a fitted model
lavPredict()
Compute estimated latent scores
lavTestLRT()
Compare (nested) lavaan models
modificationIndices()
Compute the modification indices of a model
parameterEstimates()
Parameter estimates of a latent variable model
parameterTable()
Show the table of the parameters of a fitted model
simulateData()
Simulate data starting from a lavaan model syntax
Exercise 1
Just a very easy exercise to start practicing with the lavaan sintax. The example is similar to the one above. The dataset “Exercise1.csv” contains: - N = 1100 participants - 5 variables - - gender: coded 1;2 - age: between 13 and 19 - anxiety - nevroticism - math We want to know whether anxiety and nevroticism have an effect on math achievement. Additionally, is there any effect of demographics variables?
lavaan 0.6-19 ended normally after 1 iteration
Estimator ML
Optimization method NLMINB
Number of model parameters 5
Number of observations 723
Model Test User Model:
Test statistic 0.000
Degrees of freedom 0
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
[...]
WE JUST ANALYZED NEUROTICISM AND ANXIETY. DO WE REALLY BELIEVE THEY ARE ON THE SAME LEVEL?
Path models
Path models are just pictorial representations of theoretical relationships between variables. Representations (and simulations) should be the starting point of almost every study. Based on these representations, we can build a statistical model and estimate the theoretical paths.
This, of course, is possible in lavaan.
Path models - the previous example
For example, do we really believe that anxiety and nevroticism could stay on the same level? Isn’t there a theoretical effect of one on the other?
Path models - the previous example
For example, do we really believe that anxiety and nevroticism could stay on the same level? Isn’t there a theoretical effect of one on the other?
Path models - the previous example
Code: We can write the additional regression just by adding one line to the model (plus some additional things for indirect and total effects):
m <-"math ~ age + anxiety + nevroticism + gender anxiety ~ nevroticism"fitE1 <-sem(m, data = E1)
m <-"math ~ age + am*anxiety + nm*nevroticism + gender anxiety ~ n*nevroticism # Indirect effect ind := n*am # Total effect tot := n*am + nm"fitE1 <-sem(m, data = E1)
Path model example
We have collected data from 1100 Italian cities. Research question: to what extent do economic factors, education, and environmental policies influence the quality of life of Italian people? Are these relationships mediated by social cohesion? - gdp - education - environment - wellbeing - cohesion The model can be expressed with the following equations: \[
\begin{cases}
\text{cohesion} & = \beta_{13} X_{\text{gdp}} + \beta_{14} X_{\text{education}} + \beta_{15} X_{\text{environment}} + \zeta_1 \\
\text{wellbeing} & = \beta_{23} X_{\text{gdp}} + \beta_{24} X_{\text{education}} + \beta_{25} X_{\text{environment}} + \beta_{21} X_{\text{cohesion}} + \zeta_2
\end{cases}
\]
Script
m <-"cohesion ~ gdp + education + environmentwellbeing ~ gdp + education + environment + cohesion"fit <-sem(m, data = Dprov)
Results
There is an effect of cohesion, which is in the middle.
We want to fit this model (also calculate indirect effects):
Fit a model with indirect effects
mE2 <-"#we name (_*variable) the parameters of interestlifeSatisfaction ~ a*selfEsteem + salaryselfEsteem ~ b*parentalSupport + c*attachmentattachment ~ salaryparentalSupport ~ salary#and then define the non-model parametersindAttSelf := c*aindSuppSelf := b*a"fitE2 <-sem(mE2, data = E2) #, se = "bootstrap")
Results - estimation info
lavaan 0.6-19 ended normally after 10 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 10
Number of observations 483
Model Test User Model:
Test statistic 114.172
Degrees of freedom 4
P-value (Chi-square) 0.000
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
[...]
The indirect effects estimated with lavaan in this way are just a mere multiplication of the parameters. You can apply bootstrap procedures to obtain more robust results and errors!
*’This is the goal of a good model: reproduce, from a set of theoretical associations/effects (aka covariance matrix), the original covariance matrix.
Formally:
\[
H_0 : \hat{\Sigma}(\theta) = \Sigma
\]
where \(\Sigma\) is the true covariance matrix among model variables, \(\theta\) the parameters vector, and \(\hat{\Sigma}\) the reproduced covariance matrix.’*
Model fit
\[
H_0 : \hat{\Sigma}(\theta) = \Sigma
\]
It’s time to evaluate whether our model is capable of it. To do it, we mainly compare the two matrices and obtain different fit indices.
Fit measures
Model fit refers to the ability of a model to reproduce the original covariance matrix. Fit indexes are the tools we use to estimate how good is our model in reproducing such mutrix. Most fit indexes only work with overidentified models, but some (e.g., the residualbased ones) can work with just-identified models, as well.
\(\chi^2\) based
Alternative indexes
Incremental indexes (or relative or comparative fit indexes)
Parsimony indexes
Absolute (Standalone) indexes
Residuals
\(\chi^2\) based
A T statistics following \(\chi^2\) distribution can be obtained by multiplying the sample size with the value of the fit function. df will be equal to the amount of non-redundant information minus the number of estimated parameters (remember?). This is rarely used because of its assumptions - independent observations - the non-standardized sample cov matrix is used - N is sufficiently large - manifest endogenous variables follow a multivariate normal …and because it tends to reject \(H_0\), especially with large sample sizes.
These include: - Comparative Fit Index (CFI) - Normed Fit Index (NFI) - Tucker Lewis Index (TLI) These indices compare the user model with a baseline model (the worst model).
lavaan 0.6-19 ended normally after 10 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 10
Number of observations 483
Model Test User Model:
Test statistic 114.172
Degrees of freedom 4
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 276.610
Degrees of freedom 10
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.587
Tucker-Lewis Index (TLI) -0.033
[...]
[ = ] [ = ]
With \(\delta\) being the difference between \(\chi^2\) and \(df\) and \(\delta(\text{Saturated}) = 0\)
Parsimony indexes
These include: - Information Criteria (AIC, BIC, SABIC) - Noncentrality Parameter-Based Indexes (RNI, Mc/MFI) - RMSEA: Root Mean Square Error of Approximation These models favor parsimony and penalize complex models. Among these, the RMSEA is probably the most used index. It assess how well the model approximate the data (as opposed to assessing if it is an exact fit). It is bounded between 0.0 and 1.0, with values closer to 0 (ZERO) indicating better fit.
These include: - \(\chi^2 / *df*\) ratio - GFI: Goodness of Fit Index - AGFI: Adjusted Goodness of Fit Index - SRMR/RMR: (Standardized) Root Mean Square Residual GFI and AGFI want to be 1, while the SRMR wants to be 0!
Remember that just because your model fits the data, doesn’t allow to conclude that your model is correct nor that the data generation process follows your hypothesized paths. - With path analysis, the same fit might be achieved with opposite arrows! - Test alternative models - Errors are included in manifest variables…and in the estimates! As usual all models are false, but some are useful.
What about our model?
In Exercise 2 we fit a perfect model, all our hypotheses were confirmed, effects were significant with three stars. We were happy…
Step 5 model modification
Unfortunately, it’s time to modify the model or to accept that something was wrong in our hypotheses.
Modification indices
If you have no clue about the reason why the model doesn’t converge, statistics could help you.
While the independent variable is assumed to cause the outcome variable, it’s total effect (\(c\)) is partially/totally mediated by another intervening variable, the mediator variable.
Note that our independent causes the mediator (\(a\)), and the mediator causes the outcome (\(b\)). The independent can still cause the outcome, but the path might have changed to \(c'\).
Total effects
When we have a mediation, the effect of the independent variable (X) on the outcome is given by the sum of the direct and indirect effect.
total effect = direct effect + indirect effect OR c = c’ + ab
m <-"# RegressionsY ~ c*X + b*MM ~ a*X# Indirect effectindirect := a*b# Total effecttotal := c + a*b"# Here of course we have no c'...to avoid R errors
COMMENTS?
Assumptions
The problem with mediation analyses is that they rely on often-neglected assumptions:
No X-M Interaction: The effect of M on Y or b does not vary across levels of X.
Causal Direction: The variable M causes Y, but Y does not cause M.
Perfect Reliability in M: The reliability of M is perfect.
No Confounding: There is no variable that causes M and Y.
Other limitations
Full mediation models are saturated models and we cannot obtain fit indices
You can use BIC, AIC, or SABIC and test (and compare) different models:
no direct effect model
no effect from causal variable to the mediator
no effect from the mediator to outcome
Models with inverted arrows fit the model as well as opposite models
This is true for all equivalent models. Avoid drawing paintings and base your models on strong theoretical assumptions!
Thinking first
Comparison
cfi tli srmr rmsea
0.976 0.853 0.031 0.094
cfi tli srmr rmsea
0.976 0.853 0.031 0.094
Regressions:
Estimate Std.Err z-value P(>|z|)
math ~
anxiety 0.350 0.050 6.979 0.000
sex 0.143 0.106 1.347 0.178
height 0.034 0.048 0.699 0.484
sex ~
anxiety 0.075 0.024 3.092 0.002
height ~
sex 0.984 0.091 10.807 0.000
[...]
Regressions:
Estimate Std.Err z-value P(>|z|)
math ~
anxiety 0.350 0.050 6.979 0.000
sex 0.143 0.106 1.347 0.178
height 0.034 0.048 0.699 0.484
anxiety ~
sex 0.271 0.088 3.092 0.002
height ~
sex 0.984 0.091 10.807 0.000
[...]
Readings
About mediation analyses, you can read: - the webpages by Dave Kenny: - - The history of mediations - Mediations explained and advanced topics - Roher et al. publication…but there’s a lot of information
